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e Abstract 

7 The Darwin approximation is investigated for its possible use in simulation of electromagnetic 

8 effects in large size, high frequency capacitively coupled discharges. The approximation is utilized 

CN 

(3jT) 9 within the framework of two different fluid models which are applied to typical cases showing 

10 pronounced standing wave and skin effects. With the first model it is demonstrated that Darwin 

11 approximation is valid for treatment of such effects in the range of parameters under consideration. 

12 The second approach, a reduced nonlinear Darwin approximation-based model, shows that the 

13 electromagnetic phenomena persist in a more realistic setting. The Darwin approximation offers 

14 a simple and efficient way of carrying out electromagnetic simulations as it removes the Courant 

15 condition plaguing explicit electromagnetic algorithms and can be implemented as a straightforward 
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16 modification of electrostatic algorithms. The algorithm described here avoids iterative schemes 

17 needed for the divergence cleaning and represents a fast and efficient solver, which can be used 
is in fluid and kinetic models for self-consistent description of technical plasmas exhibiting certain 
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19 electromagnetic activity. 
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20 I. INTRODUCTION 



21 Low temperature plasmas play a crucial role in materials processing [lj. Particularly 

22 capacitive radio frequency discharges are important. Such plasmas, which were first intro- 

23 duced in the 1960s, have undergone a continuous development keeping up with the increasing 

24 requirements of the hitech industry. One trend is the attempt to increase productivity by 

25 increasing the size of the processed wafers or substrates. In microelectronics, wafer sizes of 

26 300 mm are now state of the art, the transition to 450 mm is envisioned [2]. The manu- 

27 facturing of flat panel displays or photovoltaic solar cells requires discharges in the meter 

28 range [3H5]. The size of the processing equipment obviously must meet these requirements. 

29 A second recent trend is the use of higher, often multiple frequencies of up to few hun- 

30 dred MHz believed to ensure better power coupling, higher plasma density, and thus higher 

31 productivity as well [6j Cj . 

32 These developments have implications on the physical models which can be used for 

33 realistic simulations of plasma-based manufacturing processes. One such implication is of 

34 particular importance: The simultaneous increase of driving frequency, plasma density, and 

35 discharge size puts the plasmas into a regime where the commonly used "electrostatic ap- 

36 proximation" of Maxwell's equations no longer suffices. That means that V x E = (and 

37 thus E = — V$) is not justified anymore. (In this regime the skin depth and/or wavelength 

38 of the surface waves propagating in the sheath regions may become comparable or smaller 

39 than the reactor dimension). Recent experimental studies have revealed that strong plasma 

40 and field non-uniformities appear to be connected with electromagnetic effects [8HT3]. 

41 The theoretical aspects of electromagnetic effects in large-area, high frequency discharges 

42 were first systematically studied by Lieberman et al. [H] . They considered an axisymmetric 

43 slab model with uniform plasma density in the bulk and the matrix sheath (i.e., constant 

44 ion density and zero electron density), treating the plasma as a dielectric and assuming 

45 the sheath width to be constant in time. Also, a "natural boundary" condition for a CCP 

46 discharge with the dielectric sidewall that all the RF current flows through the discharge 

47 (so that the magnetic field at the sidewall is constant in the axial direction and no external 

48 source of electromagnetic excitation is involved) was employed. Under these conditions, they 

49 treated the problem analytically and gave criteria, based on the estimate of change in the 

50 power deposition profile produced by Ohmic heating due to these effects, for significance of 
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51 the standing wave and skin effects. In essence, a standing wave consists of surface waves, 

52 see, e.g., [T3], excited in the cavities formed in the powered (grounded) sheath, bounded by 

53 the electrode, powered (grounded) sheath interface and the sidewall; the surface waves in 

54 the grounded and powered regions are weakly coupled with each other through the bulk. 

55 The kind of skin effect meant here manifests itself in tendency of the induced electric field 

56 in an RF-driven discharge to reduce current at the axis of the discharge and increase it close 

57 to the discharge periphery [T6] . similar to the skin effect in metals conducting RF current 
sa (see, e.g., [IT]). 

59 Several detailed studies were conducted by Chabert et al. [18] who used a self-consistent 

eo transmission line model to account for the interaction of RF heating and plasma density. In 

ei recent years the number of contributions to the field has multiplied, for a review see, e.g., 

62 Ref. [19]. 

63 The analytical models mentioned above, being significantly limited by the underlying 

64 assumptions, leave many important questions unanswered. Numerical simulations should 

65 be used to describe the phenomena in more realistic situations. So far a number of nu- 
ee merical studies based on the fluid models have been conducted for the large scale CCP 

67 discharges driven at high frequencies [20H21]- All those approaches were based on the full 

68 set of Maxwell's equations, which contains electromagnetic radiation. The underlying equa- 

69 tions are hyperbolic and one has to satisfy the Courant criterion on numerical stability 

70 provided an explicit algorithm, such as the popular finite-difference time-domain (FDTD) 

71 algorithm, is used for their numerical solution. The Courant criterion forces the time step to 

72 be small and the overall simulation time very long. One can resort to implicit algorithms for 

73 the numerical solution of Maxwell's equations, however, it must be noted that the existing 

74 implicit numerical schemes are rather intricate. 

75 In this paper we propose an alternative numerical approach for studying the large scale, 

76 high frequency driven capacitively coupled discharges, which is able to describe the electro- 

77 magnetic phenomena occurring in such discharges. The approach is based on the Darwin 

78 approximation for solving the system of Maxwell's equations J2S]. The Darwin approxima- 

79 tion reduces the original system of hyperbolic equations to a set of elliptic equations by 
so neglecting the transversal part of the displacement current in Ampere's law, which is justi- 

81 fied by a scale analysis of the phenomena of interest. Since the resulting elliptic equations 

82 do not support electromagnetic radiation in vacuum, the corresponding Courant criterion 
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83 is removed and the simulation time step can be chosen as in electrostatic simulations. The 

84 required numerical scheme is a relatively simple and straightforward modification of the 

85 well-established electrostatic algorithms. It is also important that the approach proposed 

86 in this article avoids the computationally expensive iterative schemes used in the previous 

87 implementations of the Darwin approximation-based aglorithms (see Refs. [26H29] ). 

ss We demonstrate that the Darwin approximation is in very good agreement with the fully 

89 electromagnetic linear treatment of the problem described in p3] . However, realistic plasma 

90 density profiles are far from the model used in the latter work, hence next we study the 

91 standing wave and skin effect phenomena employing a nonlinear fluid model based on the 

92 Darwin approximation with more realistic plasma profiles treated in time domain. Unlike 

93 some works (see, e.g., Ref. [22]), which consider excitation of the electromagnetic fields 

94 in the discharge coming from external sources, we assume the "natural" excitation of the 

95 electromagnetic fields in the CCP discharge generated by the voltage/current that drives the 

96 discharge. Obviously, such kind of excitation is always present in the discharge. Also, our 

97 Darwin-fluid model accounts for finite electron inertia effects which is necessary in plasmas 

98 with low collisions driven at high frequencies. 

99 The paper is organized as follows: In Section II we briefly describe the assumptions of 

100 the Darwin approximation and its applicability. In Section III we describe the choice of 

101 representation for the electromagnetic fields that is particularly convenient for the numeri- 

102 cal implementation and a model geometry used in the examples studied in this paper, the 

103 boundary conditions are also given in this section. In Section IV we use a linear fluid model 

104 solving for the electromagnetic fields in frequency domain both with the Darwin approxi- 

105 mation and the full system of Maxwell's equations for two example cases with pronounced 
we standing wave and skin effects, respectively, to corroborate that the Darwin approximation 
107 gives an accurate description of the phenomena in question for the parameters of interest, 
los In Section IV we describe a nonlinear fluid model accounting for the finite electron inertia 
109 effects and also using the Darwin approximation to solve for the electromagnetic fields. We 
no demonstrate in a model simulation the standing wave effect for a case with more realis- 
iii tic plasma density profiles and sheath dynamics than in cases investigated in Section IV. 
n2 Finally, we summarize the main results of our study in Section VI. 
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113 II. THE DARWIN APPROXIMATION OF THE BOLTZMANN-MAXWELL SYS- 

114 TEM 



us In general, the electrodynamic behavior of any high frequency discharge can be described 

lie by the kinetic equation, such as Boltzmann's equation, supplemented with the system of 

H7 Maxwell's equations along with appropriate boundary conditions. Reduced fluid models 

us can be constructed from this general system by taking moments of Boltzmann's equation 

n9 and postulating closure relations. The dynamics of the particle distribution function of 

120 charged species s under influence of electromagnetic forces and collisions is determined by 

121 Boltzmann's equation, 



9f s . -. ^ , . Z s q s ^ g\ _ , df s 

-m +v - vfs + ^7{ E + vxB J- VJs= ^t 



(i) 



col 



122 In turn, the electromagnetic field is described by Maxwell's equations 



^VxB=J +£o f, (2) 

V • B = 0, (4) 
e V • E = p. (5) 

123 Direct attempts at numerical solution of ([l])-([5]) usually require very long simulation time. 

124 This is connected to the fact that this system supports electromagnetic waves in vacuum, 

125 which must be resolved in explicit algorithms in order to avoid numerical instabilities. An- 

126 other approach is to use an implicit algorithm for the numerical solution, where unresolved 

127 electromagnetic modes of no interest are damped numerically. However, such algorithms are 

128 cumbersome to implement and often resort to artificial numerical constructions (see, e.g., 

129 [IS]). 

130 A conceptually different alternative is to remove the stiffest time scale in electromagnetic 

131 simulation by reduction of the original system of equations using the assumption about the 

132 time and spatial scales of the expected phenomena. Namely, we use the fact that for the far 

133 most CCP discharges of interest the electrode size is small compared to the wavelength of the 

134 electromagnetic wave in vacuum corresponding to the highest frequency driving harmonics. 

135 Following the logic of [30], one can recast ^ and ^ into a normalized form where 



136 e = L/ct with L and t the typical spatial and time scale of the problem, 

VxB = j + e—, (6) 

Vx^- £ f. (7) 

137 To the zeroth order in e (J7J) yields V x E = 0, so that to this order the electric field is purely 

138 longitudinal (irrotational); we will denote this part of the electric field El- To the same 

139 order the displacement current drops out of To describe the electromagnetic effects one 

140 needs to include the next order in e, whereby 

Vx5=J + £ f, (8) 
dB 

V x Et = -e— , (9) 

141 with Et the transversal (solenoidal) part of the electric field, V • Et = 0. (|8]) and ^ 

142 constitute the Darwin approximation. 

143 One can see that this model eliminates electromagnetic waves, but keeps an important 

144 part of the physics, namely electromagnetic effects in the relatively low frequency phenomena 

145 (when the corresponding wavelength in vacuum is larger than the system size). Note that 
we the corresponding frequencies are called (very) high frequencies in the standard terminology 

147 of the CCP discharges. It is important to note that by dropping the transversal part of the 

148 displacement current the continuity equation remains satisfied. 

149 The relation of the Darwin model to the full system of Maxwell's equations can also be 

150 illustrated by the cold unmagnetized plasma dispersion relation, 

c 2 fc 2 + - u, 2 = 0. (10) 

loo + v y ' 

151 Fig. 1 displays the solutions k = k(u) of this relation for a typical parameter combination 

152 (w pc = 2tt x 1 GHz, v/ujpe = 0.01). Three distinct regimes can be discerned. Slow phenomena 



153 (oj < v) are governed by the collisional skin effect k — (1 — i) / \/2^Juj/v , with A sc f = 

154 c/ojp e the collisionless skin-depth. In the intermediate frequency range [y < uj < u pe ), 

155 the collisionless skin effect prevails, k = iX~^ { . Finally, for large frequencies (u > u pc ) this 

156 dispersion relation describes propagation of the electromagnetic waves having k = u/c. 

157 In contrast, the dispersion relation of the Darwin model is 

c 2 k 2 + ^ = 0. 11 
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158 The Darwin model then covers the middle ground between the electrostatic approximation 

159 and a fully electromagnetic treatment. 

wo Historically, Darwin's original derivation |25j was designed to describe a set of charged 

lei particles in free space with velocities small compared to the speed of light. He started from 

162 the Lagrangian description of the iV-body problem with the interaction incorporated via the 

163 retarded Lienard-Wichert potentials. A formal expansion in the smallness parameter v/c, 

164 carried up to the second order, resulted in a simplified Lagrangian with instantaneous (not 

165 retarded) potentials. In the same expansion order, relativistic corrections to the equations 
we of motion appeared. For a recent analysis of that approach, see [31 J. 
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167 III. PROBLEM SETUP 



we For the sake of concreteness, we will consider the same problem setup for studying the 

169 standing wave and skin effect phenomena in CCP discharges as was suggested by Lieberman 

170 et al. [Ej, i.e., a CCP discharge having cylindrical geometry symmetric in the azimuthal 

171 direction bounded by a dielectric sidewall (see Fig. 2). 

172 As all the numerical models described in this article will use the same field equations 

173 stemming from the Darwin approximation of Maxwell's equation, we will describe them in 

174 this section. We saw in the previous section that the electric field in the Darwin approxi- 

175 mation is separated into two parts, a longitudinal part, which is curl-free, and a transversal 

176 part, which is divergence-free. Hence, we choose to represent the electric field in the following 

177 way, 

E = -V0 + VxA T , (12) 

178 where the first (second) term expresses the longitudinal (transversal) part of the electric 

179 field. The longitudinal potential 4> is governed by Poisson's equation, 

V 2 = --, (13) 

wo with p the charge density. To find an equation for the transversal vector-potential At we 
i8i first take the curl of Ampere's law in the Darwin approximation, which results in 

vif^oVxf (14) 
where we used Q. Then, one can exploit Faraday's law using the chosen representation for 



182 



J J 5 



the transversal field, 



-» dB 

V x V x A T = — k-. (15) 
at 

84 Following [14] . we consider only TM modes with E = e r E r + e z E z and B = (IqBq, and 

85 At = egAxe, so that V • B = and V • At = hold automatically due to the symmetry of 

86 the problem. 

87 In all the following models we study the electromagnetic response to the "natural exci- 
se tation" of a CCP discharge, prescribing a harmonically varying constant magnetic field at 
89 the sidewall, which is in agreement with the previously made assumption that all the RF 
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190 current, which drives the discharge, flows through it, and no external excitation source is 

191 involved. Correspondingly, boundary conditions for Bg is d z B e = at the electrodes (z = —I 

192 and z = I), which is needed in order for the tangential electric field to vanish at the elec- 

193 trodes, and reads Bg = hqI / '(2ttR) at the sidewall (r = R) (provided no RF current leaves 

194 the discharge through the dielectric sidewall). Boundary conditions for dtBg are obtained 

195 by taking time derivative of the boundary conditions for Bg. Further, boundary conditions 
we for A T g are d z A Te = at the electrodes, so that the tangential component of the transversal 

197 electric field vanishes, and A T g = at the sidewall. The latter boundary condition for the 

198 A Te at the sidewall can be obtained as follows. From the condition that the radial compo- 

199 nent of the transversal electric field vanishes at the dielectric sidewall one can deduce that 

200 d z Axg = at r = R. Integrating this equation with respect to z at r = R, matching the 

201 boundary conditions for Are at the electrodes and setting the common constant to zero, one 

202 obtains the boundary condition sought. 

203 As the boundary conditions for the longitudinal potential will be slightly different in the 

204 two models studied in this work, we will specify the boundary conditions for each particular 

205 model separately. 
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206 IV. LINEAR MODEL PROBLEM WITH FIXED SHEATH TREATED IN FRE- 

207 QUENCY DOMAIN 



208 In this section we will demonstrate that the Darwin approximation is a valid approach 

209 for investigation of the electromagnetic phenomena in CCPs by comparing the full elec- 

210 tromagnetic solution to solution obtained with the Darwin approximation for the problem 

211 investigated by Lieberman et al. |14J. The model plasma studied there has uniform ion 

212 density and a stepwise density for the electron component, such that electrons have the 

213 same density as ions in the bulk and zero density in the sheath, the sheath boundaries are 

214 considered to be stationary. The plasma is modeled as a dielectric, change of the electric 

215 field in the sheaths due to motion of the bulk electron plasma component is modeled through 

216 the surface charge created at the interfaces between the bulk plasma and the sheath regions. 

217 This approach is valid for small oscillations of the bulk plasma (see [15j). a more general 

218 case of moving sheath boundary is considered in the next section. 



Then, using the representation of the electric field given in (12) (which holds in a general 
case as well), one can obtain the following system of equations, 

V x —V x B = x —V x A T (16) 
V x V x At = -iujB (17) 



221 provided all quantities are proportional to e luJt . In (16) and (17), €l = 1 in the sheath and 

222 6l — 1 — ujp e /u{u! — iv) in the bulk. This system of equations describes both fully EM case 

223 (if €t = €l) and the Darwin approximation (if e-r = in the sheath and €t = —(Ji e /oj(oj—ip) 

224 in the bulk; this expression can be easily obtained in the usual way by using the equation 

225 of motion and Ampere's law in the Darwin approximation). 



226 Once one obtains At and B from solution of (16) and (17), the transversal electric field 

227 can be obtained from Et = V x At and the longitudinal electric field from Ampere's law, 

228 which provides E L = —i&fujeiS x B — ^tI^iS x A t . 

229 Results in Figs. 3 and 4 show the power deposition profiles for two selected cases with 

230 similar parameters to the corresponding cases in [H]. In all cases d = 2 cm s = 0.4 cm, 

231 and v = 10 7 s -1 . The radial profiles shown in Figs. 3 and 4 are normalized values of 

232 P cap = J d \E z \ 2 dz and P ind = \E r \ 2 dz [H]. 

233 calculated with the full system of Maxwell's equations and the Darwin approximation. 
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234 Not only do both models demonstrate same behavior qualitatively, but also quantitative 

235 difference lies within a numerical error for the parameters chosen. Consequently, from this 

236 moment on we will suppose that the Darwin approximation is a valid approach for description 

237 of the standing wave and skin effects. 
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238 V. NONLINEAR MODEL PROBLEM WITH MOVING SHEATH TREATED IN 

239 TIME DOMAIN 



240 Although the model problem used [H] and in the previous section gives a certain insight 

241 into the standing wave and skin effects, it is still quite far from realistic CCP plasmas. The 

242 assumptions of stationary sheath boundary and stepwise electron density profile are the 

243 main drawbacks of the model. Most of the previous analytical and numerical treatments of 

244 the standing wave and skin effects in the literature also assumed stationary sheath boundary 

245 with the electrostatic sheath model (e.g., [22]), which makes the corresponding numerical 

246 algorithm faster, but strips the model of adequate description of possible nonlinear phe- 

247 nomena. In [21] the authors studied the problem entirely in the time domain, however, the 

248 commonly adopted drift-diffusion approximation also used in that work is clearly not appli- 

249 cable to the usual range of parameters in low-pressure CCP discharges, when the frequency 

250 of electron collisional momentum transfer is smaller than the driving frequency. In partic- 

251 ular, the highly collisional regime studied in [21] renders direct comparison with results of 

252 [14] impossible, as the latter work was made under the low-collisionality ansatz. 

253 To fully resolve the nonlinear discharge dynamics we propose another model, where the 

254 sheath boundary is allowed to move, the plasma density profiles have a more realistic shape, 

255 and the effects of finite electron mass are retained, which makes the model capable of more 

256 general treatment of the electromagnetic phenomena under question in low-pressure CCP 

257 discharges than the model considered in the previous section. The electromagnetic fields in 

258 the sheath region are calculated self-consistently with the plasma dynamics and are treated 

259 in the same way as the rest of the electromagnetic fields in the entire discharge. To keep 

260 the model simple yet comprising all the effects of interest, we adopt a fluid description for 

261 the electron plasma component, whereas ions are assumed to be immobile. In this way we 

262 choose to focus on the coupled dynamics of electrons and the electromagnetic fields in the 

263 entire discharge, omitting the issue of discharge sustainment. Analogously, we substitute 

264 the plasma energy balance equation with the assumption of electron isothermality. The 

265 geometry of the model is the same as considered in the previous section. 

266 The evolution of electrons is governed by the momentum equation, 

Ov T V?7> 6 / -> -A 

+ (v e ■ V)v e = -V0 + V x A T + v e x B - v m v e , (18) 

at m e rtp m P V / 
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267 and the continuity equation, 



dt 



-V • (n P Vp 



(19) 



268 The boundary conditions for (18) is that the radial component of the electron velocity, v e 



vanishes at the sidewall, which comes from the assumption that all the electrons are reflected 
back into the discharge there, and that the axial component of the electron velocity at the 
driven (grounded) electrodes is v £:Z = ±l/A\J%kBT e /iim e , respectively, which corresponds 



to the kinetically limited flux. The boundary conditions for (19) are, in accordance with 



the boundary conditions for (18), that the radial flux to the sidewall vanishes, and the 
axial flux to the electrodes is kinetically limited, so that n e v eiZ = ±n e /4^8kBT e /irm e at 
the driven (grounded) electrodes, respectively. We take the initial ion density profile to be 
of the Gaussian form and the initial electron density profile to be of the same form as the 
ion density profile. The numerical discharge simulation gradually develops electron-depleted 
sheath regions (see Fig. 5). 

The fields are calculated self-consistently with evolution of the electron plasma component 



using (13) to (15)). One can distinguish two different sets of boundary conditions depending 



on the way the discharge is driven, either with a voltage source or a current source. Whereas 



boundary conditions for ( 15 ) are same for both cases (see the discussion at the end of section 



283 III), the boundary conditions for (13) and (14) are different for each case 



284 A. Voltage-driven discharge 

285 In this case the potential at the driven electrode is prescribed, <p(z = —I) = V(t). This 

286 gives Dirichlet boundary condition at the driven electrode, the boundary conditions at the 

287 sidewall and the grounded electrode being same in both cases (Neumann boundary condition 

288 at the sidewall, d r <p = and at the grounded electrode Dirichlet boundary condition = 0. 

289 (See also the discussion in Section III). 

290 The time derivative of the total axial current through the discharge, which enters the 



29i sidewall boundary condition for (14), must be calculated under the circumstances by inte- 



292 grating the plasma current density time derivative, —edt(n e v e:Z ) using (18) and (19) over the 

293 surface of any discharge cross- section where the plasma current is dominant, for example, 

294 in the midplane. 
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B. Current-driven discharge 



296 The potential at the driven electrode is calculate in this case from Kirchhoff's law of 

297 current continuity at the driven electrode, where the sum of the electron current and the 

298 displacement current in the sheath must match the current supplied by the source. We 

299 seek solution to in the form <fi = (0o + <P{t))&v + &p, where V 2 $ v = 0, where 0o is the 

300 stationary, and <fi the harmonically changing parts of the potential at the driven electrode, 

301 respectively. The boundary conditions are: $ v = 1 at the driven electrode, $ v = at the 

302 grounded electrode, and <9 r $ v = at the sidewall; $ p satisfies V 2 $ p = —p/e with boundary 

303 conditions $ p = at both electrodes, and <9 r $ p = at the sidewall. From the current 

304 balance at the driven electrode one has — 2ir J^(en e v e + eod t d z (j))rdr = I(t). Substituing the 

305 chosen form for <p into this equation, one obtains an equation for the potential at the driven 

306 electrode 

f R 

I(t) + 2tt / (en e v £tZ + e d t d z ® p ) rdr 



I- ^ <*> 



27re f (d z $ v )rdr 
Jo 



307 The boundary condition for (14) at the sidewall is particularly simple in the case of 

308 current-driven discharge, as one can directly use the value of the time derivative of the total 

309 current supplied by the source. 



Finally, we want to briefly discuss some details of numerical implementation of (13)) 



311 (14), (15), (18), and (19). We use a leapfrog-like scheme to integrate (18) and (19) in time, 

312 for which we use grids staggered in time for velocity and density, the density and field values 

313 being taken at integral time levels and velocity at the half- time levels. Then, the order of 
3u the whole numerical scheme is the following: 



315 1. (18) is solved treating the nonlinear second term on the left-hand-side semi-implicitly 

316 thereby linearizing this term with respect to the density taken at the new time level. 

317 The velocity on the right hand side is time-centralized by taking the average between 

318 the new and old time levels, so that it is also semi-implicit. It is worth mentioning 

319 that during the stage of the sheath formation an artificial viscosity term is helpful 

320 in getting rid of numerical instabilities arising at the plasma-sheath interface. Once 

321 smooth sheaths are formed, the viscosity term can be switched off. 
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2. 



The density on the right-hand side of (19) is treated semi-implicitly and is time- 
centralized using the same technique as in the previous procedure. 



3. (13) is solved with the boundary conditions depending on the way the discharge is 



driven according to the techniques described above. 



326 4. The source term in (14) should be calculated using the expression d t j e = —e(v e d t n e + 



n e d t v e ) with d t v e and d t n e calculated with help of (18) and (19). An attempt at direct 
evaluation of d t j e will lead to numerical instability because numerical discretization 
of a time derivative introduces numerical dispersion and finite field propagation time, 
whereas conceptually the Darwin approximation expects instantaneous field change in 



the whole domain (see, e.g., [32]). When calculating the source term for (14), it is also 
useful to write the term /x V x e 2 n e /m e V x A T as — X^ cl d t B + fx e 2 n e /m e Vn e x V x A T , 



where A sc f = c/uj pe is the collisionless skin depth and we used (15). The first term in 
this expression describes the collisionless skin effect and can be transferred to the left 
hand side to be treated implicitly for better numerical stability. 



336 5. Finally, (15) is solved using the boundary conditions discussed in Section III. 



337 One can see that because of the chosen representation for the transversal field as Et = V x 

338 At the cumbersome and computationally expensive iterative scheme needed for making sure 

339 that the calculated Et is divergence-free used in the previous implementations of the Darwin 

340 codes (see, e.g., [26H29] ) is avoided. This is because Et calculated by taking the curl of At 

341 guarantees that Et is divergence-free automatically. To account for the electromagnetic 

342 effects in 2D geometry of the problem we described, one needs to solve only two additional 

343 elliptic equations using the same field solver that is used for solving Poisson's equation. It is 

344 a small price for removing the very expensive Courant criterion connected with propagation 

345 of the light wave, though. 

346 In Fig. 5 we show the radial profile of relative power deposition due to the Ohmic heat- 

347 ing resulting from a voltage-driven discharge simulation showing a significant standing wave 

348 effect for a case with relatively low collisionality {y m = 1 x 10 7 s _1 ). Parameters are sim- 

349 ilar to the parameters for the case shown on the top part of Fig. 3. The inductive power 

350 deposition is negligible and is not shown. The change of radial power deposition profile is 

351 more pronounced for the more realistic Gaussian-like axial electron density profile adopted 
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352 for the simulation in contrast to the stepwise electron density profile taken for the problem 

353 in Section IV. For comparison, we also present the radial profile of relative power deposition 

354 obtained with a purely electrostatic simulation. As one can see, the electromagnetic compo- 

355 nent is an essential part of the standing wave effect for the parameters in question, despite 

356 the fact that surface waves propagating in sheath and generating standing waves may also 

357 exist in the electrostatic description (see, e.g., [33J). 
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358 VI. CONCLUSIONS AND OUTLOOK 



359 In this paper we have demonstrated that the Darwin approximation is adequate for de- 

360 scription of the standing wave and skin effects present in the high frequency CCP discharges 

361 by comparing relevant electromagnetic fields calculated using the approximation and solu- 

362 tion of the full set of Maxwell's equations for two different typical cases exhibiting such 

363 effects in a model problem. We considered the "natural" excitation of such effects coming 

364 from the current driving the discharge needed for its very existence. The Darwin approxima- 

365 tion reduces originally hyperbolic Maxwell's equations, by dropping the transversal current 

366 in Ampere's law, to a set of elliptic equations which do not support light waves propagating 
36? in vacuum. This has very significant consequence for the numerical algorithms based on 

368 that approach that the corresponding Courant criterion is also removed, enabling same time 

369 step size as in electrostatic algorithms. We have shown that the Darwin approximation can 

370 be straightforwardly implemented in fluid models, requiring only a slight modification of 

371 a possibly existing electrostatic code with the eventual great computational pay-off. The 

372 choice of the representation for the transversal field proposed by us helps to avoid the com- 

373 putationally demanding iterative schemes needed for its divergence cleaning used in other 

374 implementations. Using a reduced nonlinear fluid code self-consistently evolving plasma and 

375 electromagnetic fields in time domain, on the example of a case exhibiting standing wave 

376 effect we have shown that the electromagnetic effects not only are present but can be even 

377 more significant for more realistic profiles and sheath dynamics. 

378 The proper description of low-pressure CCP discharges ultimately requires a kinetic tool 

379 as kinetic effects are essential in such discharges. The use of the Darwin approximation is, 

380 of course, not limited to fluid models. The sheath dynamics, the heating mechanisms other 

381 than the Ohmic heating, the nonlocal transport phenomena all demand kinetic treatment. 

382 The Darwin numerical algorithm for the solution of the electromagnetic fields designed in 

383 this work can be also used in a Particle-in-Cell code. The reduced fluid model developed in 

384 this work can, in turn, serve as benchmark for such a code. 
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FIG. 1. The dispersion relation obtained with Darwin approximation closely follows the full elec- 
tromagnetic dispersion relation as long as u < u pe . Solid curves denote the corresponding real and 
dashed ones imagn 
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FIG. 2. Sketch of the problem geometry studied in both fluid models. The cylindrical CCP 
discharge under scrutiny has symmetric electrodes, a sidewall made of dielectric and a plasma. 
The boundary at the sheath-plasma interface is stationary in the model described in Section IV, 
whereas in the model of Section V it is allowed to move under the influence of electromagnetic 
fields. The electromagnetic fields under interest consist of the TM modes in such a plasma-filled 
cavity and have nonzero E z ,E r components for the electric, and Bq component for the magnetic 
fields, respectively. 
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FIG. 3. Comparison of solutions obtained with Darwin approximation (DA) and full system of 
Maxwell's equations (EM). The first case, which exhibits a strong standing wave effect, has / = 
40, 7 MHz and n = 10 15 uT 3 . 
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FIG. 4. The second case shows a strong skin effect and has / = 13.56 MHz and n 
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FIG. 5. Development of the sheath regions with depleted electrons in a typical simulation 
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FIG. 6. Relative power deposition radial profile (Ohmic heating) for a case with / = 40.7 MHz 
and n eymax = 1 x 10 15 m -3 
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